
      subroutine block(nvtx, ijkvtx, nsp, itdim, iwid, jwid, kwid, c1x, c2x, &
         c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, &
         c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvolaxis, vol, vvol, clite, dt&
         , qom, qdnv, cntr, bxv, byv, bzv, beta1, beta2, beta3, a11, a12, a13, &
         a21, a22, a23, a31, a32, a33)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: nsp
      integer , intent(in) :: itdim
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      real(double) , intent(in) :: clite
      real(double) , intent(in) :: dt
      real(double) , intent(in) :: cntr
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vvolaxis(*)
      real(double) , intent(in) :: vol(*)
      real(double)  :: vvol(*)
      real(double) , intent(in) :: qom(*)
      real(double) , intent(in) :: qdnv(itdim,*)
      real(double) , intent(in) :: bxv(*)
      real(double) , intent(in) :: byv(*)
      real(double) , intent(in) :: bzv(*)
      real(double) , intent(inout) :: beta1(*)
      real(double) , intent(inout) :: beta2(*)
      real(double) , intent(inout) :: beta3(*)
      real(double) , intent(out) :: a11(*)
      real(double) , intent(out) :: a12(*)
      real(double) , intent(out) :: a13(*)
      real(double) , intent(out) :: a21(*)
      real(double) , intent(out) :: a22(*)
      real(double) , intent(out) :: a23(*)
      real(double) , intent(out) :: a31(*)
      real(double) , intent(out) :: a32(*)
      real(double) , intent(out) :: a33(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, is
      real(double) :: dtsq, cdtsq, alfa, ws, rvvol, rvimjk, rvimjmk, rvijmk, &
         rvijk, rvimjkm, rvimjmkm, rvijmkm, rvijkm, agxx, agxy, agxz, agyy, &
         agyz, agzz, vdtsq
!-----------------------------------------------
!
!
!     A ROUTINE TO CALCULATE THE ELEMENTS OF THE APPROXIMATE
!     INVERSE FOR MAXWELL'S EQUATIONS
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
         a11(ijk) = 0.0
         a12(ijk) = 0.0
         a13(ijk) = 0.0
!
         a21(ijk) = 0.0
         a22(ijk) = 0.0
         a23(ijk) = 0.0
!
         a31(ijk) = 0.0
         a32(ijk) = 0.0
         a33(ijk) = 0.0
!
         beta1(ijk) = 0.0
         beta2(ijk) = 0.0
         beta3(ijk) = 0.0
!
      end do
!
      dtsq = 0.5*cntr*dt**2
      cdtsq = 2.0*clite**2*cntr*dtsq
!test      CDTSQ=0.0
      do is = 1, nsp
!
         alfa = 0.5*qom(is)*dt/clite
!
!dir$ ivdep
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            ws = qom(is)*qdnv(ijk,is)/(1. + alfa**2*(bxv(ijk)**2+byv(ijk)**2+&
               bzv(ijk)**2))
!
            beta1(ijk) = beta1(ijk) + ws
            beta2(ijk) = beta2(ijk) + ws*alfa
            beta3(ijk) = beta3(ijk) + ws*alfa**2
!
         end do
      end do
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
!      rvvol=vvolaxis(ijk)/(vvol(ijk)+1.e-20)
         rvvol = 1.0
!
         rvimjk = 1./(vol(ijk-iwid)+1.E-20)
         rvimjmk = 1./(vol(ijk-iwid-jwid)+1.E-20)
         rvijmk = 1./(vol(ijk-jwid)+1.E-20)
         rvijk = 1./(vol(ijk)+1.E-20)
!
         rvimjkm = 1./(vol(ijk-iwid-kwid)+1.E-20)
         rvimjmkm = 1./(vol(ijk-iwid-jwid-kwid)+1.E-20)
         rvijmkm = 1./(vol(ijk-jwid-kwid)+1.E-20)
         rvijkm = 1./(vol(ijk-kwid)+1.E-20)
!
         agxx = -(c1x(ijk-iwid)*c1x(ijk-iwid)*rvimjk+c2x(ijk-iwid-jwid)*c2x(ijk&
            -iwid-jwid)*rvimjmk+c3x(ijk-jwid)*c3x(ijk-jwid)*rvijmk+c4x(ijk)*c4x&
            (ijk)*rvijk+c5x(ijk-iwid-kwid)*c5x(ijk-iwid-kwid)*rvimjkm+c6x(ijk-&
            iwid-jwid-kwid)*c6x(ijk-iwid-jwid-kwid)*rvimjmkm+c7x(ijk-jwid-kwid)&
            *c7x(ijk-jwid-kwid)*rvijmkm+c8x(ijk-kwid)*c8x(ijk-kwid)*rvijkm)*&
            rvvol
!
         agxy = -(c1x(ijk-iwid)*c1y(ijk-iwid)*rvimjk+c2x(ijk-iwid-jwid)*c2y(ijk&
            -iwid-jwid)*rvimjmk+c3x(ijk-jwid)*c3y(ijk-jwid)*rvijmk+c4x(ijk)*c4y&
            (ijk)*rvijk+c5x(ijk-iwid-kwid)*c5y(ijk-iwid-kwid)*rvimjkm+c6x(ijk-&
            iwid-jwid-kwid)*c6y(ijk-iwid-jwid-kwid)*rvimjmkm+c7x(ijk-jwid-kwid)&
            *c7y(ijk-jwid-kwid)*rvijmkm+c8x(ijk-kwid)*c8y(ijk-kwid)*rvijkm)*&
            rvvol
!
         agxz = -(c1x(ijk-iwid)*c1z(ijk-iwid)*rvimjk+c2x(ijk-iwid-jwid)*c2z(ijk&
            -iwid-jwid)*rvimjmk+c3x(ijk-jwid)*c3z(ijk-jwid)*rvijmk+c4x(ijk)*c4z&
            (ijk)*rvijk+c5x(ijk-iwid-kwid)*c5z(ijk-iwid-kwid)*rvimjkm+c6x(ijk-&
            iwid-jwid-kwid)*c6z(ijk-iwid-jwid-kwid)*rvimjmkm+c7x(ijk-jwid-kwid)&
            *c7z(ijk-jwid-kwid)*rvijmkm+c8x(ijk-kwid)*c8z(ijk-kwid)*rvijkm)*&
            rvvol
!
         agyy = -(c1y(ijk-iwid)*c1y(ijk-iwid)*rvimjk+c2y(ijk-iwid-jwid)*c2y(ijk&
            -iwid-jwid)*rvimjmk+c3y(ijk-jwid)*c3y(ijk-jwid)*rvijmk+c4y(ijk)*c4y&
            (ijk)*rvijk+c5y(ijk-iwid-kwid)*c5y(ijk-iwid-kwid)*rvimjkm+c6y(ijk-&
            iwid-jwid-kwid)*c6y(ijk-iwid-jwid-kwid)*rvimjmkm+c7y(ijk-jwid-kwid)&
            *c7y(ijk-jwid-kwid)*rvijmkm+c8y(ijk-kwid)*c8y(ijk-kwid)*rvijkm)*&
            rvvol
!
         agyz = -(c1y(ijk-iwid)*c1z(ijk-iwid)*rvimjk+c2y(ijk-iwid-jwid)*c2z(ijk&
            -iwid-jwid)*rvimjmk+c3y(ijk-jwid)*c3z(ijk-jwid)*rvijmk+c4y(ijk)*c4z&
            (ijk)*rvijk+c5y(ijk-iwid-kwid)*c5z(ijk-iwid-kwid)*rvimjkm+c6y(ijk-&
            iwid-jwid-kwid)*c6z(ijk-iwid-jwid-kwid)*rvimjmkm+c7y(ijk-jwid-kwid)&
            *c7z(ijk-jwid-kwid)*rvijmkm+c8y(ijk-kwid)*c8z(ijk-kwid)*rvijkm)*&
            rvvol
!
         agzz = -(c1z(ijk-iwid)*c1z(ijk-iwid)*rvimjk+c2z(ijk-iwid-jwid)*c2z(ijk&
            -iwid-jwid)*rvimjmk+c3z(ijk-jwid)*c3z(ijk-jwid)*rvijmk+c4z(ijk)*c4z&
            (ijk)*rvijk+c5z(ijk-iwid-kwid)*c5z(ijk-iwid-kwid)*rvimjkm+c6z(ijk-&
            iwid-jwid-kwid)*c6z(ijk-iwid-jwid-kwid)*rvimjmkm+c7z(ijk-jwid-kwid)&
            *c7z(ijk-jwid-kwid)*rvijmkm+c8z(ijk-kwid)*c8z(ijk-kwid)*rvijkm)*&
            rvvol
!
         vdtsq = vvolaxis(ijk)*dtsq
!test      vdtsq=0.
!
         a11(ijk) = vvolaxis(ijk) - (agxx + agyy + agzz)*cdtsq + vdtsq*(beta1(&
            ijk)+beta3(ijk)*bxv(ijk)**2)
!
         a12(ijk) = vdtsq*(beta2(ijk)*bzv(ijk)+beta3(ijk)*byv(ijk)*bxv(ijk))
!
         a13(ijk) = vdtsq*((-beta2(ijk)*byv(ijk))+beta3(ijk)*bzv(ijk)*bxv(ijk))
!
         a21(ijk) = vdtsq*((-beta2(ijk)*bzv(ijk))+beta3(ijk)*bxv(ijk)*byv(ijk))
!
         a22(ijk) = vvolaxis(ijk) - (agxx + agyy + agzz)*cdtsq + vdtsq*(beta1(&
            ijk)+beta3(ijk)*byv(ijk)**2)
!
         a23(ijk) = vdtsq*(beta2(ijk)*bxv(ijk)+beta3(ijk)*bzv(ijk)*byv(ijk))
!
         a31(ijk) = vdtsq*(beta2(ijk)*byv(ijk)+beta3(ijk)*bxv(ijk)*bzv(ijk))
!
         a32(ijk) = vdtsq*((-beta2(ijk)*bxv(ijk))+beta3(ijk)*byv(ijk)*bzv(ijk))
!
         a33(ijk) = vvolaxis(ijk) - (agxx + agyy + agzz)*cdtsq + vdtsq*(beta1(&
            ijk)+beta3(ijk)*bzv(ijk)**2)
!
!
!
      end do
!
      return
      end subroutine block


      subroutine curlc(ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
         c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
         c3z, c4z, c5z, c6z, c7z, c8z, vol, vx, vy, vz, curlcx, curlcy, curlcz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkcell(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vol(*)
      real(double) , intent(in) :: vx(*)
      real(double) , intent(in) :: vy(*)
      real(double) , intent(in) :: vz(*)
      real(double) , intent(out) :: curlcx(*)
      real(double) , intent(out) :: curlcy(*)
      real(double) , intent(out) :: curlcz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: rvol
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, ncells
!
         ijk = ijkcell(n)
!
         rvol = 1./vol(ijk)
!
         curlcx(ijk) = (c1y(ijk)*vz(ijk+iwid)+c2y(ijk)*vz(ijk+iwid+jwid)+c3y(&
            ijk)*vz(ijk+jwid)+c4y(ijk)*vz(ijk)+c5y(ijk)*vz(ijk+iwid+kwid)+c6y(&
            ijk)*vz(ijk+iwid+jwid+kwid)+c7y(ijk)*vz(ijk+jwid+kwid)+c8y(ijk)*vz(&
            ijk+kwid)-c1z(ijk)*vy(ijk+iwid)-c2z(ijk)*vy(ijk+iwid+jwid)-c3z(ijk)&
            *vy(ijk+jwid)-c4z(ijk)*vy(ijk)-c5z(ijk)*vy(ijk+iwid+kwid)-c6z(ijk)*&
            vy(ijk+iwid+jwid+kwid)-c7z(ijk)*vy(ijk+jwid+kwid)-c8z(ijk)*vy(ijk+&
            kwid))*rvol
!
         curlcy(ijk) = (c1z(ijk)*vx(ijk+iwid)+c2z(ijk)*vx(ijk+iwid+jwid)+c3z(&
            ijk)*vx(ijk+jwid)+c4z(ijk)*vx(ijk)+c5z(ijk)*vx(ijk+iwid+kwid)+c6z(&
            ijk)*vx(ijk+iwid+jwid+kwid)+c7z(ijk)*vx(ijk+jwid+kwid)+c8z(ijk)*vx(&
            ijk+kwid)-c1x(ijk)*vz(ijk+iwid)-c2x(ijk)*vz(ijk+iwid+jwid)-c3x(ijk)&
            *vz(ijk+jwid)-c4x(ijk)*vz(ijk)-c5x(ijk)*vz(ijk+iwid+kwid)-c6x(ijk)*&
            vz(ijk+iwid+jwid+kwid)-c7x(ijk)*vz(ijk+jwid+kwid)-c8x(ijk)*vz(ijk+&
            kwid))*rvol
!
         curlcz(ijk) = (c1x(ijk)*vy(ijk+iwid)+c2x(ijk)*vy(ijk+iwid+jwid)+c3x(&
            ijk)*vy(ijk+jwid)+c4x(ijk)*vy(ijk)+c5x(ijk)*vy(ijk+iwid+kwid)+c6x(&
            ijk)*vy(ijk+iwid+jwid+kwid)+c7x(ijk)*vy(ijk+jwid+kwid)+c8x(ijk)*vy(&
            ijk+kwid)-c1y(ijk)*vx(ijk+iwid)-c2y(ijk)*vx(ijk+iwid+jwid)-c3y(ijk)&
            *vx(ijk+jwid)-c4y(ijk)*vx(ijk)-c5y(ijk)*vx(ijk+iwid+kwid)-c6y(ijk)*&
            vx(ijk+iwid+jwid+kwid)-c7y(ijk)*vx(ijk+jwid+kwid)-c8y(ijk)*vx(ijk+&
            kwid))*rvol
!
      end do
!
      return
      end subroutine curlc


      subroutine curlv(nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, vx, vy, vz, curlvx, curlvy, curlvz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vx(*)
      real(double) , intent(in) :: vy(*)
      real(double) , intent(in) :: vz(*)
      real(double) , intent(out) :: curlvx(*)
      real(double) , intent(out) :: curlvy(*)
      real(double) , intent(out) :: curlvz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
!
         curlvx(ijk) = -(c1y(ijk-iwid)*vz(ijk-iwid)+c2y(ijk-iwid-jwid)*vz(ijk-&
            iwid-jwid)+c3y(ijk-jwid)*vz(ijk-jwid)+c4y(ijk)*vz(ijk)+c5y(ijk-iwid&
            -kwid)*vz(ijk-iwid-kwid)+c6y(ijk-iwid-jwid-kwid)*vz(ijk-iwid-jwid-&
            kwid)+c7y(ijk-jwid-kwid)*vz(ijk-jwid-kwid)+c8y(ijk-kwid)*vz(ijk-&
            kwid)-c1z(ijk-iwid)*vy(ijk-iwid)-c2z(ijk-iwid-jwid)*vy(ijk-iwid-&
            jwid)-c3z(ijk-jwid)*vy(ijk-jwid)-c4z(ijk)*vy(ijk)-c5z(ijk-iwid-kwid&
            )*vy(ijk-iwid-kwid)-c6z(ijk-iwid-jwid-kwid)*vy(ijk-iwid-jwid-kwid)-&
            c7z(ijk-jwid-kwid)*vy(ijk-jwid-kwid)-c8z(ijk-kwid)*vy(ijk-kwid))
!
         curlvy(ijk) = -(c1z(ijk-iwid)*vx(ijk-iwid)+c2z(ijk-iwid-jwid)*vx(ijk-&
            iwid-jwid)+c3z(ijk-jwid)*vx(ijk-jwid)+c4z(ijk)*vx(ijk)+c5z(ijk-iwid&
            -kwid)*vx(ijk-iwid-kwid)+c6z(ijk-iwid-jwid-kwid)*vx(ijk-iwid-jwid-&
            kwid)+c7z(ijk-jwid-kwid)*vx(ijk-jwid-kwid)+c8z(ijk-kwid)*vx(ijk-&
            kwid)-c1x(ijk-iwid)*vz(ijk-iwid)-c2x(ijk-iwid-jwid)*vz(ijk-iwid-&
            jwid)-c3x(ijk-jwid)*vz(ijk-jwid)-c4x(ijk)*vz(ijk)-c5x(ijk-iwid-kwid&
            )*vz(ijk-iwid-kwid)-c6x(ijk-iwid-jwid-kwid)*vz(ijk-iwid-jwid-kwid)-&
            c7x(ijk-jwid-kwid)*vz(ijk-jwid-kwid)-c8x(ijk-kwid)*vz(ijk-kwid))
!
         curlvz(ijk) = -(c1x(ijk-iwid)*vy(ijk-iwid)+c2x(ijk-iwid-jwid)*vy(ijk-&
            iwid-jwid)+c3x(ijk-jwid)*vy(ijk-jwid)+c4x(ijk)*vy(ijk)+c5x(ijk-iwid&
            -kwid)*vy(ijk-iwid-kwid)+c6x(ijk-iwid-jwid-kwid)*vy(ijk-iwid-jwid-&
            kwid)+c7x(ijk-jwid-kwid)*vy(ijk-jwid-kwid)+c8x(ijk-kwid)*vy(ijk-&
            kwid)-c1y(ijk-iwid)*vx(ijk-iwid)-c2y(ijk-iwid-jwid)*vx(ijk-iwid-&
            jwid)-c3y(ijk-jwid)*vx(ijk-jwid)-c4y(ijk)*vx(ijk)-c5y(ijk-iwid-kwid&
            )*vx(ijk-iwid-kwid)-c6y(ijk-iwid-jwid-kwid)*vx(ijk-iwid-jwid-kwid)-&
            c7y(ijk-jwid-kwid)*vx(ijk-jwid-kwid)-c8y(ijk-kwid)*vx(ijk-kwid))
!
      end do
!
      return
      end subroutine curlv


      subroutine diagonal(ncells, ijkcell, tiny, nstart, iwid, jwid, kwid, c1x&
         , c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y&
         , c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvol, vol, diag)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: nstart
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      real(double) , intent(in) :: tiny
      integer , intent(in) :: ijkcell(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vvol(*)
      real(double)  :: vol(*)
      real(double) , intent(out) :: diag(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
!-----------------------------------------------
!
!dir$ ivdep
      do n = nstart, ncells
!
         ijk = ijkcell(n)
!
         diag(ijk) = ((c1x(ijk)**2+c1y(ijk)**2+c1z(ijk)**2)/vvol(ijk+iwid)+(c2x&
            (ijk)**2+c2y(ijk)**2+c2z(ijk)**2)/vvol(ijk+iwid+jwid)+(c3x(ijk)**2+&
            c3y(ijk)**2+c3z(ijk)**2)/vvol(ijk+jwid)+(c4x(ijk)**2+c4y(ijk)**2+&
            c4z(ijk)**2)/vvol(ijk)+(c5x(ijk)**2+c5y(ijk)**2+c5z(ijk)**2)/vvol(&
            ijk+iwid+kwid)+(c6x(ijk)**2+c6y(ijk)**2+c6z(ijk)**2)/vvol(ijk+iwid+&
            jwid+kwid)+(c7x(ijk)**2+c7y(ijk)**2+c7z(ijk)**2)/vvol(ijk+jwid+kwid&
            )+(c8x(ijk)**2+c8y(ijk)**2+c8z(ijk)**2)/vvol(ijk+kwid)) - 1.*tiny
!
      end do
!
!dir$ ivdep
      do n = 1, nstart - 1
!
         ijk = ijkcell(n)
!
         diag(ijk) = ((c1x(ijk)**2+c1y(ijk)**2+c1z(ijk)**2)/vvol(ijk+iwid)+(c2x&
            (ijk)**2+c2y(ijk)**2+c2z(ijk)**2)/vvol(ijk+iwid+jwid)+(c3x(ijk)**2+&
            c3y(ijk)**2+c3z(ijk)**2)/vvol(ijk+jwid)+(c4x(ijk)**2+c4y(ijk)**2+&
            c4z(ijk)**2)/vvol(ijk)+(c5x(ijk)**2+c5y(ijk)**2+c5z(ijk)**2)/vvol(&
            ijk+iwid+kwid)+(c6x(ijk)**2+c6y(ijk)**2+c6z(ijk)**2)/vvol(ijk+iwid+&
            jwid+kwid)+(c7x(ijk)**2+c7y(ijk)**2+c7z(ijk)**2)/vvol(ijk+jwid+kwid&
            )+(c8x(ijk)**2+c8y(ijk)**2+c8z(ijk)**2)/vvol(ijk+kwid))
!
      end do
!
      return
      end subroutine diagonal


      subroutine divc(ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
         c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
         c3z, c4z, c5z, c6z, c7z, c8z, vol, ex, ey, ez, dive)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkcell(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vol(*)
      real(double) , intent(in) :: ex(*)
      real(double) , intent(in) :: ey(*)
      real(double) , intent(in) :: ez(*)
      real(double) , intent(out) :: dive(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: divex, divey, divez
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, ncells
!
         ijk = ijkcell(n)
!
         divex = c1x(ijk)*ex(ijk+iwid) + (c2x(ijk)*ex(ijk+iwid+jwid)+(c3x(ijk)*&
            ex(ijk+jwid)+(c4x(ijk)*ex(ijk)+(c5x(ijk)*ex(ijk+iwid+kwid)+(c6x(ijk&
            )*ex(ijk+iwid+jwid+kwid)+(c7x(ijk)*ex(ijk+jwid+kwid)+c8x(ijk)*ex(&
            ijk+kwid)))))))
!
         divey = c1y(ijk)*ey(ijk+iwid) + (c2y(ijk)*ey(ijk+iwid+jwid)+(c3y(ijk)*&
            ey(ijk+jwid)+(c4y(ijk)*ey(ijk)+(c5y(ijk)*ey(ijk+iwid+kwid)+(c6y(ijk&
            )*ey(ijk+iwid+jwid+kwid)+(c7y(ijk)*ey(ijk+jwid+kwid)+c8y(ijk)*ey(&
            ijk+kwid)))))))
!
         divez = c1z(ijk)*ez(ijk+iwid) + (c2z(ijk)*ez(ijk+iwid+jwid)+(c3z(ijk)*&
            ez(ijk+jwid)+(c4z(ijk)*ez(ijk)+(c5z(ijk)*ez(ijk+iwid+kwid)+(c6z(ijk&
            )*ez(ijk+iwid+jwid+kwid)+(c7z(ijk)*ez(ijk+jwid+kwid)+c8z(ijk)*ez(&
            ijk+kwid)))))))
!
         dive(ijk) = (divex + divey + divez)/vol(ijk)
!
      end do
!
      return
      end subroutine divc


      subroutine divv(ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, &
         c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, &
         c3z, c4z, c5z, c6z, c7z, c8z, vx, vy, vz, diverge)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkcell(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vx(*)
      real(double) , intent(in) :: vy(*)
      real(double) , intent(in) :: vz(*)
      real(double) , intent(out) :: diverge(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: divergex, divergey, divergez
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, ncells
!
         ijk = ijkcell(n)
!
!
!
         divergex = -(c1x(ijk-iwid)*vx(ijk-iwid)+c2x(ijk-iwid-jwid)*vx(ijk-iwid&
            -jwid)+c3x(ijk-jwid)*vx(ijk-jwid)+c4x(ijk)*vx(ijk)+c5x(ijk-iwid-&
            kwid)*vx(ijk-iwid-kwid)+c6x(ijk-iwid-jwid-kwid)*vx(ijk-iwid-jwid-&
            kwid)+c7x(ijk-jwid-kwid)*vx(ijk-jwid-kwid)+c8x(ijk-kwid)*vx(ijk-&
            kwid))
!
         divergey = -(c1y(ijk-iwid)*vy(ijk-iwid)+c2y(ijk-iwid-jwid)*vy(ijk-iwid&
            -jwid)+c3y(ijk-jwid)*vy(ijk-jwid)+c4y(ijk)*vy(ijk)+c5y(ijk-iwid-&
            kwid)*vy(ijk-iwid-kwid)+c6y(ijk-iwid-jwid-kwid)*vy(ijk-iwid-jwid-&
            kwid)+c7y(ijk-jwid-kwid)*vy(ijk-jwid-kwid)+c8y(ijk-kwid)*vy(ijk-&
            kwid))
!
         divergez = -(c1z(ijk-iwid)*vz(ijk-iwid)+c2z(ijk-iwid-jwid)*vz(ijk-iwid&
            -jwid)+c3z(ijk-jwid)*vz(ijk-jwid)+c4z(ijk)*vz(ijk)+c5z(ijk-iwid-&
            kwid)*vz(ijk-iwid-kwid)+c6z(ijk-iwid-jwid-kwid)*vz(ijk-iwid-jwid-&
            kwid)+c7z(ijk-jwid-kwid)*vz(ijk-jwid-kwid)+c8z(ijk-kwid)*vz(ijk-&
            kwid))
!
         diverge(ijk) = divergex + divergey + divergez
!
      end do
!
      return
      end subroutine divv


      subroutine divpi(nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, pixx, pixy, pixz, piyy, piyz, pizz, divpix&
         , divpiy, divpiz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: pixx(*)
      real(double) , intent(in) :: pixy(*)
      real(double) , intent(in) :: pixz(*)
      real(double) , intent(in) :: piyy(*)
      real(double) , intent(in) :: piyz(*)
      real(double) , intent(in) :: pizz(*)
      real(double) , intent(out) :: divpix(*)
      real(double) , intent(out) :: divpiy(*)
      real(double) , intent(out) :: divpiz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: divxx, divyx, divzx, divxy, divyy, divzy, divxz, divyz, &
         divzz
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
!
         divxx = -(c1x(ijk-iwid)*pixx(ijk-iwid)+(c2x(ijk-iwid-jwid)*pixx(ijk-&
            iwid-jwid)+(c3x(ijk-jwid)*pixx(ijk-jwid)+(c4x(ijk)*pixx(ijk)+(c5x(&
            ijk-iwid-kwid)*pixx(ijk-iwid-kwid)+(c6x(ijk-iwid-jwid-kwid)*pixx(&
            ijk-iwid-jwid-kwid)+(c7x(ijk-jwid-kwid)*pixx(ijk-jwid-kwid)+c8x(ijk&
            -kwid)*pixx(ijk-kwid))))))))
!
         divyx = -(c1y(ijk-iwid)*pixy(ijk-iwid)+(c2y(ijk-iwid-jwid)*pixy(ijk-&
            iwid-jwid)+(c3y(ijk-jwid)*pixy(ijk-jwid)+(c4y(ijk)*pixy(ijk)+(c5y(&
            ijk-iwid-kwid)*pixy(ijk-iwid-kwid)+(c6y(ijk-iwid-jwid-kwid)*pixy(&
            ijk-iwid-jwid-kwid)+(c7y(ijk-jwid-kwid)*pixy(ijk-jwid-kwid)+c8y(ijk&
            -kwid)*pixy(ijk-kwid))))))))
!
         divzx = -(c1z(ijk-iwid)*pixz(ijk-iwid)+(c2z(ijk-iwid-jwid)*pixz(ijk-&
            iwid-jwid)+(c3z(ijk-jwid)*pixz(ijk-jwid)+(c4z(ijk)*pixz(ijk)+(c5z(&
            ijk-iwid-kwid)*pixz(ijk-iwid-kwid)+(c6z(ijk-iwid-jwid-kwid)*pixz(&
            ijk-iwid-jwid-kwid)+(c7z(ijk-jwid-kwid)*pixz(ijk-jwid-kwid)+c8z(ijk&
            -kwid)*pixz(ijk-kwid))))))))
!
         divpix(ijk) = divxx + divyx + divzx
!
      end do
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
         divxy = -(c1x(ijk-iwid)*pixy(ijk-iwid)+(c2x(ijk-iwid-jwid)*pixy(ijk-&
            iwid-jwid)+(c3x(ijk-jwid)*pixy(ijk-jwid)+(c4x(ijk)*pixy(ijk)+(c5x(&
            ijk-iwid-kwid)*pixy(ijk-iwid-kwid)+(c6x(ijk-iwid-jwid-kwid)*pixy(&
            ijk-iwid-jwid-kwid)+(c7x(ijk-jwid-kwid)*pixy(ijk-jwid-kwid)+c8x(ijk&
            -kwid)*pixy(ijk-kwid))))))))
!
         divyy = -(c1y(ijk-iwid)*piyy(ijk-iwid)+(c2y(ijk-iwid-jwid)*piyy(ijk-&
            iwid-jwid)+(c3y(ijk-jwid)*piyy(ijk-jwid)+(c4y(ijk)*piyy(ijk)+(c5y(&
            ijk-iwid-kwid)*piyy(ijk-iwid-kwid)+(c6y(ijk-iwid-jwid-kwid)*piyy(&
            ijk-iwid-jwid-kwid)+(c7y(ijk-jwid-kwid)*piyy(ijk-jwid-kwid)+c8y(ijk&
            -kwid)*piyy(ijk-kwid))))))))
!
         divzy = -(c1z(ijk-iwid)*piyz(ijk-iwid)+(c2z(ijk-iwid-jwid)*piyz(ijk-&
            iwid-jwid)+(c3z(ijk-jwid)*piyz(ijk-jwid)+(c4z(ijk)*piyz(ijk)+(c5z(&
            ijk-iwid-kwid)*piyz(ijk-iwid-kwid)+(c6z(ijk-iwid-jwid-kwid)*piyz(&
            ijk-iwid-jwid-kwid)+(c7z(ijk-jwid-kwid)*piyz(ijk-jwid-kwid)+c8z(ijk&
            -kwid)*piyz(ijk-kwid))))))))
!
         divpiy(ijk) = divxy + divyy + divzy
!
      end do
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
         divxz = -(c1x(ijk-iwid)*pixz(ijk-iwid)+(c2x(ijk-iwid-jwid)*pixz(ijk-&
            iwid-jwid)+(c3x(ijk-jwid)*pixz(ijk-jwid)+(c4x(ijk)*pixz(ijk)+(c5x(&
            ijk-iwid-kwid)*pixz(ijk-iwid-kwid)+(c6x(ijk-iwid-jwid-kwid)*pixz(&
            ijk-iwid-jwid-kwid)+(c7x(ijk-jwid-kwid)*pixz(ijk-jwid-kwid)+c8x(ijk&
            -kwid)*pixz(ijk-kwid))))))))
!
         divyz = -(c1y(ijk-iwid)*piyz(ijk-iwid)+(c2y(ijk-iwid-jwid)*piyz(ijk-&
            iwid-jwid)+(c3y(ijk-jwid)*piyz(ijk-jwid)+(c4y(ijk)*piyz(ijk)+(c5y(&
            ijk-iwid-kwid)*piyz(ijk-iwid-kwid)+(c6y(ijk-iwid-jwid-kwid)*piyz(&
            ijk-iwid-jwid-kwid)+(c7y(ijk-jwid-kwid)*piyz(ijk-jwid-kwid)+c8y(ijk&
            -kwid)*piyz(ijk-kwid))))))))
!
         divzz = -(c1z(ijk-iwid)*pizz(ijk-iwid)+(c2z(ijk-iwid-jwid)*pizz(ijk-&
            iwid-jwid)+(c3z(ijk-jwid)*pizz(ijk-jwid)+(c4z(ijk)*pizz(ijk)+(c5z(&
            ijk-iwid-kwid)*pizz(ijk-iwid-kwid)+(c6z(ijk-iwid-jwid-kwid)*pizz(&
            ijk-iwid-jwid-kwid)+(c7z(ijk-jwid-kwid)*pizz(ijk-jwid-kwid)+c8z(ijk&
            -kwid)*pizz(ijk-kwid))))))))
!
         divpiz(ijk) = divxz + divyz + divzz
!
      end do
!
      return
      end subroutine divpi


      subroutine gradc(nvtx, ijkvtx, kbp1, iwid, jwid, kwid, c1x, c2x, c3x, c4x&
         , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z&
         , c3z, c4z, c5z, c6z, c7z, c8z, vol, psi, gradcx, gradcy, gradcz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer  :: kbp1
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vol(*)
      real(double) , intent(in) :: psi(*)
      real(double) , intent(out) :: gradcx(*)
      real(double) , intent(out) :: gradcy(*)
      real(double) , intent(out) :: gradcz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: rvol
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
!
         rvol = 1./vol(ijk)
!
         gradcx(ijk) = (c1x(ijk)*psi(ijk+iwid)+c2x(ijk)*psi(ijk+iwid+jwid)+c3x(&
            ijk)*psi(ijk+jwid)+c4x(ijk)*psi(ijk)+c5x(ijk)*psi(ijk+iwid+kwid)+&
            c6x(ijk)*psi(ijk+iwid+jwid+kwid)+c7x(ijk)*psi(ijk+jwid+kwid)+c8x(&
            ijk)*psi(ijk+kwid))*rvol
!
         gradcy(ijk) = (c1y(ijk)*psi(ijk+iwid)+c2y(ijk)*psi(ijk+iwid+jwid)+c3y(&
            ijk)*psi(ijk+jwid)+c4y(ijk)*psi(ijk)+c5y(ijk)*psi(ijk+iwid+kwid)+&
            c6y(ijk)*psi(ijk+iwid+jwid+kwid)+c7y(ijk)*psi(ijk+jwid+kwid)+c8y(&
            ijk)*psi(ijk+kwid))*rvol
!
         gradcz(ijk) = (c1z(ijk)*psi(ijk+iwid)+c2z(ijk)*psi(ijk+iwid+jwid)+c3z(&
            ijk)*psi(ijk+jwid)+c4z(ijk)*psi(ijk)+c5z(ijk)*psi(ijk+iwid+kwid)+&
            c6z(ijk)*psi(ijk+iwid+jwid+kwid)+c7z(ijk)*psi(ijk+jwid+kwid)+c8z(&
            ijk)*psi(ijk+kwid))*rvol
!
      end do
!
      return
      end subroutine gradc



      subroutine dnull(ncells, ijkcell, iwid, jwid, kwid, nsmooth, f, dive, vol&
         )
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: nsmooth
      integer , intent(in) :: ijkcell(*)
      real(double) , intent(inout) :: f(*)
      real(double) , intent(inout) :: dive(*)
      real(double)  :: vol(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: iter, n, ijk
      real(double) :: tiny
!-----------------------------------------------
!
      do iter = 1, nsmooth
!
         tiny = 0.01
!
!dir$ ivdep
         do n = 1, ncells
!
            ijk = ijkcell(n)
!
!     eps=tiny/vol(ijk)**2
!
            dive(ijk) = 8.*f(ijk) - 4.*(f(ijk+iwid)+f(ijk-iwid)+f(ijk+jwid)+f(&
               ijk-jwid)+f(ijk+kwid)+f(ijk-kwid)) + 2.*(f(ijk+iwid+jwid)+f(ijk+&
               iwid-jwid)+f(ijk-iwid+jwid)+f(ijk-iwid-jwid)+f(ijk+iwid+kwid)+f(&
               ijk+iwid-kwid)+f(ijk-iwid+kwid)+f(ijk-iwid-kwid)+f(ijk+kwid+jwid&
               )+f(ijk+kwid-jwid)+f(ijk-kwid+jwid)+f(ijk-kwid-jwid)) - (f(ijk+&
               iwid+jwid+kwid)+f(ijk+iwid+jwid-kwid)+f(ijk+iwid-jwid+kwid)+f(&
               ijk+iwid-jwid-kwid)+f(ijk-iwid+jwid+kwid)+f(ijk-iwid+jwid-kwid)+&
               f(ijk-iwid-jwid+kwid)+f(ijk-iwid-jwid-kwid))
!
         end do
!dir$ ivdep
         do n = 1, ncells
!
            ijk = ijkcell(n)
!
            f(ijk) = (1. - tiny)*f(ijk) + dive(ijk)*tiny
!
         end do
!
      end do
!
      return
      end subroutine dnull




      subroutine gnull(nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, f, gradxf, gradyf, gradzf, vvol)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double) , intent(in) :: f(*)
      real(double) , intent(inout) :: gradxf(*)
      real(double) , intent(inout) :: gradyf(*)
      real(double) , intent(inout) :: gradzf(*)
      real(double) , intent(in) :: vvol(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: tiny, eps, gtiny
!-----------------------------------------------
!
      tiny = 0.1
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
!      eps=tiny*vvol(ijk)**(2./3.)*(mod(n,2)*2-1)
         eps = tiny*vvol(ijk)**(2./3.)
!
         gtiny = eps*(f(ijk-iwid)-f(ijk-iwid-jwid)+f(ijk-jwid)-f(ijk)-f(ijk-&
            iwid-kwid)+f(ijk-iwid-jwid-kwid)-f(ijk-jwid-kwid)+f(ijk-kwid))/vvol&
            (ijk)
         gradxf(ijk) = gradxf(ijk) + gtiny
         gradyf(ijk) = gradyf(ijk) + gtiny
         gradzf(ijk) = gradzf(ijk) + gtiny
!
      end do
!
      return
      end subroutine gnull




      subroutine gradf(nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x&
         , c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z&
         , c4z, c5z, c6z, c7z, c8z, f, gradxf, gradyf, gradzf)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: f(*)
      real(double) , intent(out) :: gradxf(*)
      real(double) , intent(out) :: gradyf(*)
      real(double) , intent(out) :: gradzf(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
         gradxf(ijk) = -(c1x(ijk-iwid)*f(ijk-iwid)+(c2x(ijk-iwid-jwid)*f(ijk-&
            iwid-jwid)+(c3x(ijk-jwid)*f(ijk-jwid)+(c4x(ijk)*f(ijk)+(c5x(ijk-&
            iwid-kwid)*f(ijk-iwid-kwid)+(c6x(ijk-iwid-jwid-kwid)*f(ijk-iwid-&
            jwid-kwid)+(c7x(ijk-jwid-kwid)*f(ijk-jwid-kwid)+c8x(ijk-kwid)*f(ijk&
            -kwid))))))))
!
         gradyf(ijk) = -(c1y(ijk-iwid)*f(ijk-iwid)+(c2y(ijk-iwid-jwid)*f(ijk-&
            iwid-jwid)+(c3y(ijk-jwid)*f(ijk-jwid)+(c4y(ijk)*f(ijk)+(c5y(ijk-&
            iwid-kwid)*f(ijk-iwid-kwid)+(c6y(ijk-iwid-jwid-kwid)*f(ijk-iwid-&
            jwid-kwid)+(c7y(ijk-jwid-kwid)*f(ijk-jwid-kwid)+c8y(ijk-kwid)*f(ijk&
            -kwid))))))))
!
         gradzf(ijk) = -(c1z(ijk-iwid)*f(ijk-iwid)+(c2z(ijk-iwid-jwid)*f(ijk-&
            iwid-jwid)+(c3z(ijk-jwid)*f(ijk-jwid)+(c4z(ijk)*f(ijk)+(c5z(ijk-&
            iwid-kwid)*f(ijk-iwid-kwid)+(c6z(ijk-iwid-jwid-kwid)*f(ijk-iwid-&
            jwid-kwid)+(c7z(ijk-jwid-kwid)*f(ijk-jwid-kwid)+c8z(ijk-kwid)*f(ijk&
            -kwid))))))))
!
      end do
!
      return
      end subroutine gradf


      subroutine graddiv(ncells, ijkcell, iwid, jwid, kwid, ibp1, jbp1, kbp1, &
         cdlt, sdlt, periodicx, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, ex, ey, ez, divetil, gdivx, gdivy, gdivz)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      real(double)  :: cdlt
      real(double)  :: sdlt
      logical  :: periodicx
      integer  :: ijkcell(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double)  :: divetil(*)
      real(double)  :: gdivx(*)
      real(double)  :: gdivy(*)
      real(double)  :: gdivz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
!-----------------------------------------------
!
!
!       ************************************************************************
!
!     CALCULATE GRAD(DIV(E))
!
      call list (1, ibp1, 1, jbp1, 1, kbp1, iwid, jwid, kwid, ncells, ijkcell)
!
      call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, vol, ex, ey, ez, divetil)
!
      call list (2, ibp1, 2, jbp1, 2, kbp1, iwid, jwid, kwid, ncells, ijkcell)
!
      call gradf (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, divetil, gdivx, gdivy, gdivz)
!
      call bcperv (ibp1 + 1, jbp1 + 1, kbp1 + 1, iwid, jwid, kwid, sdlt, cdlt, &
         gdivx, gdivy, gdivz, periodicx)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gdivx, gdivy, gdivz)
!
!
!     END CALCULATION OF GRAD(DIV(E))
!
!     ********************************************************************
!
      return
      end subroutine graddiv


      subroutine list(i1, i2, j1, j2, k1, k2, iwid, jwid, kwid, nlist, ijklist)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: i1
      integer , intent(in) :: i2
      integer , intent(in) :: j1
      integer , intent(in) :: j2
      integer , intent(in) :: k1
      integer , intent(in) :: k2
      integer , intent(in) :: iwid
      integer , intent(in) :: jwid
      integer , intent(in) :: kwid
      integer , intent(out) :: nlist
      integer , intent(out) :: ijklist(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: k, j, i
!-----------------------------------------------
!
      nlist = 0
!
      do k = k1, k2
         do j = j1, j2
            do i = i1, i2
!
               nlist = nlist + 1
!
               ijklist(nlist) = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid
!
            end do
         end do
      end do
!
      return
      end subroutine list


      subroutine ludecomp(nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32&
         , a33)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(inout) :: a11(*)
      real(double) , intent(inout) :: a12(*)
      real(double) , intent(inout) :: a13(*)
      real(double) , intent(inout) :: a21(*)
      real(double) , intent(inout) :: a22(*)
      real(double) , intent(inout) :: a23(*)
      real(double) , intent(inout) :: a31(*)
      real(double) , intent(inout) :: a32(*)
      real(double) , intent(inout) :: a33(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: b32, b33
!-----------------------------------------------
!
!dir$ ivdep
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
         a11(ijk) = 1./a11(ijk)
!
         a21(ijk) = a21(ijk)*a11(ijk)
         a22(ijk) = a22(ijk) - a21(ijk)*a12(ijk)
         a23(ijk) = a23(ijk) - a21(ijk)*a13(ijk)
!
         a31(ijk) = a31(ijk)*a11(ijk)
         b32 = a32(ijk) - a31(ijk)*a12(ijk)
         b33 = a33(ijk) - a31(ijk)*a13(ijk)
!
         a22(ijk) = 1./a22(ijk)
!
         a32(ijk) = b32*a22(ijk)
         a33(ijk) = b33 - a32(ijk)*a23(ijk)
!
         a33(ijk) = 1./a33(ijk)
!
         a23(ijk) = a23(ijk)*a22(ijk)
         a12(ijk) = a12(ijk)*a11(ijk)
         a13(ijk) = a13(ijk)*a11(ijk)
!
      end do
!
      return
      end subroutine ludecomp


      subroutine mudote(nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt&
         , transpos, clite, ex, ey, ez, exmu, eymu, ezmu)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: nsp
      integer , intent(in) :: itdim
      real(double) , intent(in) :: dt
      real(double) , intent(in) :: transpos
      real(double) , intent(in) :: clite
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: qdnv(itdim,*)
      real(double) , intent(in) :: qom(*)
      real(double) , intent(in) :: bxv(*)
      real(double) , intent(in) :: byv(*)
      real(double) , intent(in) :: bzv(*)
      real(double) , intent(in) :: ex(*)
      real(double) , intent(in) :: ey(*)
      real(double) , intent(in) :: ez(*)
      real(double) , intent(inout) :: exmu(*)
      real(double) , intent(inout) :: eymu(*)
      real(double) , intent(inout) :: ezmu(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, is
      real(double) :: beta, omcx, omcy, omcz, edotb, denom
!-----------------------------------------------
!
!dir$ ivdep
!
      exmu(ijkvtx(:nvtx)) = 0.0
      eymu(ijkvtx(:nvtx)) = 0.0
      ezmu(ijkvtx(:nvtx)) = 0.0
!
!
      do is = 1, nsp
!
         beta = 0.5*qom(is)*dt*transpos/clite
!
!dir$ ivdep
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            omcx = beta*bxv(ijk)
            omcy = beta*byv(ijk)
            omcz = beta*bzv(ijk)
!
            edotb = ex(ijk)*omcx + (ey(ijk)*omcy+ez(ijk)*omcz)
!
            denom = qom(is)*qdnv(ijk,is)/(1. + (omcx**2 + omcy**2 + omcz**2))
!
            exmu(ijk) = exmu(ijk) + (ex(ijk)+(ey(ijk)*omcz-ez(ijk)*omcy+edotb*&
               omcx))*denom
!
            eymu(ijk) = eymu(ijk) + (ey(ijk)+(ez(ijk)*omcx-ex(ijk)*omcz+edotb*&
               omcy))*denom
!
            ezmu(ijk) = ezmu(ijk) + (ez(ijk)+(ex(ijk)*omcy-ey(ijk)*omcx+edotb*&
               omcz))*denom
!
         end do
      end do
!
      return
      end subroutine mudote
!
      subroutine sustensor(nvtx,ijkvtx,nsp, itdim, qdnv, qom, bxv, byv, bzv, dt&
         , transpos, clite, cntr, susxx, susxy, susxz,suszx, suszy, suszz)
! ... to calculate the susceptibility tensor: epsilon = I + mu
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: nsp
      integer , intent(in) :: itdim
      real(double) , intent(in) :: dt
      real(double) , intent(in) :: transpos
      real(double) , intent(in) :: clite
      real(double) , intent(in) :: cntr
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: qdnv(itdim,*)
      real(double) , intent(in) :: qom(*)
      real(double) , intent(in) :: bxv(*)
      real(double) , intent(in) :: byv(*)
      real(double) , intent(in) :: bzv(*)
      real(double) , intent(inout) :: susxx(*)
      real(double) , intent(inout) :: susxy(*)
      real(double) , intent(inout) :: susxz(*)
!      real(double) , intent(inout) :: susyx(*)
!      real(double) , intent(inout) :: susyy(*)
!      real(double) , intent(inout) :: susyz(*)
      real(double) , intent(inout) :: suszx(*)
      real(double) , intent(inout) :: suszy(*)
      real(double) , intent(inout) :: suszz(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, is
      real(double) :: beta, omcx, omcy, omcz, denom
!-----------------------------------------------
!
!dir$ ivdep
!
      susxx(ijkvtx(:nvtx)) = 1.0
      susxy(ijkvtx(:nvtx)) = 0.0
      susxz(ijkvtx(:nvtx)) = 0.0
!      susyx(ijkvtx(:nvtx)) = 0.0
!      susyy(ijkvtx(:nvtx)) = 1.0
!      susyz(ijkvtx(:nvtx)) = 0.0
      suszx(ijkvtx(:nvtx)) = 0.0
      suszy(ijkvtx(:nvtx)) = 0.0
      suszz(ijkvtx(:nvtx)) = 1.0
!
!
      do is = 1, nsp
!
         beta = 0.5*qom(is)*dt*transpos/clite
!
!dir$ ivdep
         do n = 1, nvtx
!
            ijk = ijkvtx(n)
!
            omcx = beta*bxv(ijk)
            omcy = beta*byv(ijk)
            omcz = beta*bzv(ijk)
!
            denom = qom(is)*qdnv(ijk,is)/(1. + (omcx**2 + omcy**2 +&
                  omcz**2))*0.5*cntr*dt**2
!
            susxx(ijk) = susxx(ijk) + (1.d0+omcx**2)*denom
            susxy(ijk) = susxy(ijk) + (omcz+omcx*omcy)*denom
            susxz(ijk) = susxz(ijk) + (omcx*omcz-omcy)*denom
!            susyx(ijk) = susyx(ijk) + (omcx*omcy-omcz)*denom
!            susyy(ijk) = susyy(ijk) + (1.d0+omcy**2)*denom
!            susyz(ijk) = susyz(ijk) + (omcx+omcy*omcz)*denom
            suszx(ijk) = suszx(ijk) + (omcy+omcx*omcz)*denom
            suszy(ijk) = suszy(ijk) + (omcy*omcz-omcx)*denom
            suszz(ijk) = suszz(ijk) + (1.d0+omcz**2)*denom
!
         end do
      end do
!
      return
      end subroutine sustensor

!
      subroutine newerpois(ncells, ijkcell, nvtx, ijkvtx, itdim, nsp, iwid, &
         jwid, kwid, precon, transpos, tiny, ibp1, jbp1, kbp1, c1x, c2x, c3x, &
         c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, &
         c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvolaxis, itsub, itmax, &
         error, fourpi, dt, cntr, clite, wkl, wkr, wkf, wke, periodicx, qdnc, &
         qdnv, qom, bxv, byv, bzv, exmu, eymu, ezmu, ex, ey, ez, dive, residu, &
         aq, q, gradphix, gradphiy, gradphiz, phi, diag, elle, nsmooth, itertot&
         )
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtx
      integer  :: itdim
      integer  :: nsp
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: kbp1
      integer  :: itsub
      integer , intent(in) :: itmax
      integer  :: nsmooth
      integer , intent(out) :: itertot
      real(double)  :: transpos
      real(double)  :: tiny
      real(double)  :: error
      real(double)  :: fourpi
      real(double)  :: dt
      real(double)  :: cntr
      real(double)  :: clite
      real(double)  :: wkl
      real(double)  :: wkr
      real(double)  :: wkf
      real(double)  :: wke
      real(double)  :: elle
      logical  :: precon
      logical  :: periodicx
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: vvol(*)
      real(double)  :: vvolaxis(*)
      real(double)  :: qdnc(*)
      real(double)  :: qdnv(itdim,*)
      real(double)  :: qom(*)
      real(double)  :: bxv(*)
      real(double)  :: byv(*)
      real(double)  :: bzv(*)
      real(double)  :: exmu(*)
      real(double)  :: eymu(*)
      real(double)  :: ezmu(*)
      real(double)  :: ex(*)
      real(double)  :: ey(*)
      real(double)  :: ez(*)
      real(double)  :: dive(*)
      real(double)  :: residu(*)
      real(double)  :: aq(*)
      real(double)  :: q(itdim,20)
      real(double)  :: gradphix(*)
      real(double)  :: gradphiy(*)
      real(double)  :: gradphiz(*)
      real(double)  :: phi(*)
      real(double)  :: diag(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp2, jbp2, kbp2, n, ijk, it, iter
      real(double) :: phibc, bnorm, rnorm, phiavg, volavg, exavg, eyavg, ezavg&
         , eenrg, totvol, rvvol
!-----------------------------------------------
!
!
!
      phibc = 0.
      ibp2 = ibp1 + 1
      jbp2 = jbp1 + 1
      kbp2 = kbp1 + 1
!     a routine to solve poisson's equation
!
!                  -div((mu)dot(grad(phi)))=4*pi*(n-div(jtilde)*cntr*dt)
!
!     USING GMRES(itsub), itsub<=20
!
!     boundary conditions correspond to a topological torus
!     that is periodic in i and j, and with Dirichlet bc's at k=kbp2
!
!     calculate source term
!
      diag(ijkcell(:ncells)) = 0.0
!
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , diag, periodicx)
!
      call newres (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, jbp1&
         , tiny, cntr, nsp, itdim, dt, transpos, clite, qdnv, qom, bxv, byv, &
         bzv, exmu, eymu, ezmu, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, diag, gradphix, gradphiy, gradphiz, vvol, dive, qdnc, fourpi&
         , residu, elle)
!
      bnorm = 0.0
      do n = 1, ncells
         ijk = ijkcell(n)
         bnorm = bnorm + abs(residu(ijk))
      end do
      write (*, *) 'bnorm NEWERpoisson', bnorm
!
      if (bnorm < error) bnorm = 1.0
!
      if (precon) then
!
         call axisvol (ibp1, jbp1, iwid, jwid, kwid, vvol, nvtx, ijkvtx, vol, &
            vvolaxis)
!
         call preconcr (ncells, ijkcell, nsp, itdim, iwid, jwid, kwid, qdnv, &
            qom, bxv, byv, bzv, dt, transpos, clite, cntr, c1x, c2x, c3x, c4x, &
            c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, &
            c2z, c3z, c4z, c5z, c6z, c7z, c8z, vvol, diag)
!
      else
!
!dir$ ivdep
         diag(ijkcell(:ncells)) = -1.
!
      endif
!
!
      do it = 1, itmax
!
         call gmres (ncells, ijkcell, nvtx, ijkvtx, itdim, nsp, iwid, jwid, &
            kwid, precon, transpos, tiny, ibp1, jbp1, kbp1, c1x, c2x, c3x, c4x&
            , c5x, c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, &
            c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, vvol, vvolaxis, itsub, iter&
            , error, rnorm, fourpi, dt, cntr, clite, qdnc, qdnv, qom, bxv, byv&
            , bzv, exmu, eymu, ezmu, dive, residu, aq, q, bnorm, gradphix, &
            gradphiy, gradphiz, phi, diag, elle, nsmooth, wkl, wkr, wkf, wke, &
            periodicx)
!
!     Added by GL, subtracts average
!
         phiavg = 0.
         volavg = 0.
         phiavg = phiavg + sum(phi(ijkcell(:ncells))*vol(ijkcell(:ncells)))
         volavg = sum(vol(ijkcell(:ncells)))
!
!     UNCOMMENT FOR NEUMANN
!
!       phi(ijk)=phi(ijk)-phiavg
         phiavg = phiavg/volavg
         call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, &
            phibc, phi, periodicx)
!
         if (iter < itsub) go to 98
      end do
!
      write (*, *) 'gmres fails', (it - 1)*itsub + iter, rnorm
      itertot = (it - 1)*itsub + iter
!     write(*,*)ibp2*jbp2*kbp2
      phi(:ibp2*jbp2*kbp2) = 0.
      go to 99
   98 continue
      write (*, *) 'gmres converges ', (it - 1)*itsub + iter, rnorm
      itertot = (it - 1)*itsub + iter
   99 continue
      call bcphi (ibp1, jbp1, kbp1, iwid, jwid, kwid, wkl, wkr, wkf, wke, phibc&
         , phi, periodicx)
!
      call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, phi, gradphix, gradphiy, gradphiz)
!
!
!      call axisgrad(ibp1,jbp1,iwid,jwid,kwid,
!     &     gradphix,gradphiy,gradphiz)
!
!dir$ ivdep
      exavg = 0.
      eyavg = 0.
      ezavg = 0.
      eenrg = 0.
      totvol = 0.
      do n = 1, nvtx
         ijk = ijkvtx(n)
         rvvol = 1./vvol(ijk)
         ex(ijk) = -gradphix(ijk)*rvvol
         ey(ijk) = -gradphiy(ijk)*rvvol
         ez(ijk) = -gradphiz(ijk)*rvvol
         totvol = totvol + vvol(ijk)
         exavg = exavg + ex(ijk)*vvol(ijk)
         eyavg = eyavg + ey(ijk)*vvol(ijk)
         ezavg = ezavg + ez(ijk)*vvol(ijk)
         eenrg = eenrg + (ex(ijk)**2+ey(ijk)**2+ez(ijk)**2)*vvol(ijk)
      end do
!
      exavg = exavg/totvol
      eyavg = eyavg/totvol
      ezavg = ezavg/totvol
!
!
!     test solution
!
      call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, vol, ex, ey, ez, dive)
!
      bnorm = 0.0
      bnorm = dot_product(dive(ijkcell(:ncells)),vol(ijkcell(:ncells)))
      write (*, *) 'test NEWERpoisson: integral of div E: ', bnorm
!
      call newres (ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, jbp1&
         , tiny, cntr, nsp, itdim, dt, transpos, clite, qdnv, qom, bxv, byv, &
         bzv, exmu, eymu, ezmu, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, phi, gradphix, gradphiy, gradphiz, vvol, dive, qdnc, fourpi&
         , residu, elle)
!
      bnorm = 0.0
      do n = 1, ncells
         ijk = ijkcell(n)
         bnorm = bnorm + abs(residu(ijk))
      end do
      write (*, *) 'test NEWERpoisson', bnorm
!
      return
      end subroutine newerpois


!
!
      subroutine newres(ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1, &
         jbp1, tiny, cntr, nsp, itdim, dt, transpos, clite, qdnv, qom, bxv, byv&
         , bzv, exmu, eymu, ezmu, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, &
         c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, &
         c8z, vol, phi, gradphix, gradphiy, gradphiz, vvol, dive, qdnc, fourpi&
         , residu, elle)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtx
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      integer  :: nsp
      integer  :: itdim
      real(double)  :: tiny
      real(double) , intent(in) :: cntr
      real(double)  :: dt
      real(double)  :: transpos
      real(double)  :: clite
      real(double) , intent(in) :: fourpi
      real(double) , intent(in) :: elle
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: qdnv(itdim,*)
      real(double)  :: qom(*)
      real(double)  :: bxv(*)
      real(double)  :: byv(*)
      real(double)  :: bzv(*)
      real(double)  :: exmu(*)
      real(double)  :: eymu(*)
      real(double)  :: ezmu(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: phi(*)
      real(double)  :: gradphix(*)
      real(double)  :: gradphiy(*)
      real(double)  :: gradphiz(*)
      real(double) , intent(in) :: vvol(*)
      real(double)  :: dive(*)
      real(double) , intent(in) :: qdnc(*)
      real(double) , intent(out) :: residu(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: lambda, lamb1, dtsq, rvvol, tau
!-----------------------------------------------
!
!
      dtsq = dt**2
!
!     calculate residue
!
      call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, phi, gradphix, gradphiy, gradphiz)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradphix, gradphiy, gradphiz&
         )
!
      do n = 1, nvtx
         ijk = ijkvtx(n)
         rvvol = 1./vvol(ijk)
         gradphix(ijk) = -gradphix(ijk)*rvvol
         gradphiy(ijk) = -gradphiy(ijk)*rvvol
         gradphiz(ijk) = -gradphiz(ijk)*rvvol
      end do
!

      nsp = 2
      call mudote (nvtx, ijkvtx, nsp, itdim, qdnv, qom, bxv, byv, bzv, dt, &
         transpos, clite, gradphix, gradphiy, gradphiz, exmu, eymu, ezmu)
!
      do n = 1, nvtx
         ijk = ijkvtx(n)
         rvvol = 1./vvol(ijk)
         exmu(ijk) = gradphix(ijk) + 0.5*cntr*dtsq*exmu(ijk)
         eymu(ijk) = gradphiy(ijk) + 0.5*cntr*dtsq*eymu(ijk)
         ezmu(ijk) = gradphiz(ijk) + 0.5*cntr*dtsq*ezmu(ijk)
!      exmu(ijk)=gradphix(ijk)+cntr*dtsq*exmu(ijk)
!      eymu(ijk)=gradphiy(ijk)+cntr*dtsq*eymu(ijk)
!      ezmu(ijk)=gradphiz(ijk)+cntr*dtsq*ezmu(ijk)
      end do
!
!     calculate div(grad(phi))
!
      call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, vol, exmu, eymu, ezmu, dive)
!
!
!     A call to reduce null space
!
!      call dnull(ncells,ijkcell,iwid,jwid,kwid,
!     &     phi,dive,vol)
!
      tau = 0.
      lamb1 = 0. + tau
      lambda = 1.
      if (elle /= 0.) lamb1 = 1./elle
!
      do n = 1, ncells
         ijk = ijkcell(n)
!
         residu(ijk) = (fourpi*qdnc(ijk)*lambda**2-dive(ijk)-lamb1**2*phi(ijk))&
            *vol(ijk)
!
      end do
!
      return
      end subroutine newres


      subroutine ndcplr(nstart, nvtx, ijkvtx, iw, jw, kw, tiny, h, vvol, &
         delsqphi)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nstart
      integer , intent(in) :: nvtx
      integer , intent(in) :: iw
      integer , intent(in) :: jw
      integer , intent(in) :: kw
      real(double) , intent(in) :: tiny
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(inout) :: h(*)
      real(double) , intent(in) :: vvol(*)
      real(double) , intent(inout) :: delsqphi(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: r48, cplr
!-----------------------------------------------
!
!dir$ ivdep
!
      r48 = 1./48.
      do n = nstart, nvtx
!
         ijk = ijkvtx(n)
!
         cplr = ((-24.*h(ijk)*vvol(ijk))+8.*(h(ijk+iw)*vvol(ijk+iw)+h(ijk-iw)*&
            vvol(ijk-iw)+h(ijk+jw)*vvol(ijk+jw)+h(ijk-jw)*vvol(ijk-jw)+h(ijk+kw&
            )*vvol(ijk+kw)+h(ijk-kw)*vvol(ijk-kw))-2.*(h(ijk+iw+jw)*vvol(ijk+iw&
            -jw)+h(ijk-iw+jw)*vvol(ijk-iw+jw)+h(ijk-iw-jw)*vvol(ijk-iw-jw)+h(&
            ijk+iw-jw)*vvol(ijk+iw-jw)+h(ijk+jw+kw)*vvol(ijk+jw+kw)+h(ijk-jw+kw&
            )*vvol(ijk-jw+kw)+h(ijk-jw-kw)*vvol(ijk-jw-kw)+h(ijk+jw-kw)*vvol(&
            ijk+jw-kw)+h(ijk+iw+kw)*vvol(ijk+iw+kw)+h(ijk-iw+kw)*vvol(ijk-iw+kw&
            )+h(ijk+iw-kw)*vvol(ijk+iw-kw)+h(ijk-iw-kw)*vvol(ijk-iw-kw)))*r48
!
         delsqphi(ijk) = h(ijk)*vvol(ijk) + tiny*cplr
!
      end do
!
      do n = nstart, nvtx
!
         ijk = ijkvtx(n)
!
         h(ijk) = delsqphi(ijk)/vvol(ijk)
!
      end do
!
!
      return
      end subroutine ndcplr



      subroutine preconcr(ncells, ijkcell, nsp, itdim, iw, jw, kw, qdnv, qom, &
         bxv, byv, bzv, dt, transpos, clite, cntr, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, vvol, diag)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: ncells
      integer , intent(in) :: nsp
      integer , intent(in) :: itdim
      integer , intent(in) :: iw
      integer , intent(in) :: jw
      integer , intent(in) :: kw
      real(double) , intent(in) :: dt
      real(double) , intent(in) :: transpos
      real(double) , intent(in) :: clite
      real(double) , intent(in) :: cntr
      integer , intent(in) :: ijkcell(*)
      real(double) , intent(in) :: qdnv(itdim,*)
      real(double) , intent(in) :: qom(*)
      real(double) , intent(in) :: bxv(*)
      real(double) , intent(in) :: byv(*)
      real(double) , intent(in) :: bzv(*)
      real(double) , intent(in) :: c1x(*)
      real(double) , intent(in) :: c2x(*)
      real(double) , intent(in) :: c3x(*)
      real(double) , intent(in) :: c4x(*)
      real(double) , intent(in) :: c5x(*)
      real(double) , intent(in) :: c6x(*)
      real(double) , intent(in) :: c7x(*)
      real(double) , intent(in) :: c8x(*)
      real(double) , intent(in) :: c1y(*)
      real(double) , intent(in) :: c2y(*)
      real(double) , intent(in) :: c3y(*)
      real(double) , intent(in) :: c4y(*)
      real(double) , intent(in) :: c5y(*)
      real(double) , intent(in) :: c6y(*)
      real(double) , intent(in) :: c7y(*)
      real(double) , intent(in) :: c8y(*)
      real(double) , intent(in) :: c1z(*)
      real(double) , intent(in) :: c2z(*)
      real(double) , intent(in) :: c3z(*)
      real(double) , intent(in) :: c4z(*)
      real(double) , intent(in) :: c5z(*)
      real(double) , intent(in) :: c6z(*)
      real(double) , intent(in) :: c7z(*)
      real(double) , intent(in) :: c8z(*)
      real(double) , intent(in) :: vvol(*)
      real(double) , intent(out) :: diag(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk, is
      real(double) :: cdotb1, cdotb2, cdotb3, cdotb4, cdotb5, cdotb6, cdotb7, &
         dtcsq, wscx1, wscy1, wscz1, wscx2, wscy2, wscz2, wscx3, wscy3, wscz3, &
         wscx4, wscy4, wscz4, wscx5, wscy5, wscz5, wscx6, wscy6, wscz6, wscx7, &
         wscy7, wscz7, wscx8, wscy8, wscz8, beta, cdotb, denom, divwscx, &
         divwscy, divwscz
!-----------------------------------------------
!
!
      dtcsq = 0.5*cntr*dt**2
!
      do n = 1, ncells
!
         ijk = ijkcell(n)
!
         wscx1 = -c1x(ijk)/vvol(ijk+iw)
         wscy1 = -c1y(ijk)/vvol(ijk+iw)
         wscz1 = -c1z(ijk)/vvol(ijk+iw)
!/vvol(ijk)
         wscx2 = -c2x(ijk)/vvol(ijk+iw+jw)
         wscy2 = -c2y(ijk)/vvol(ijk+iw+jw)
         wscz2 = -c2z(ijk)/vvol(ijk+iw+jw)
!/vvol(ijk)
         wscx3 = -c3x(ijk)/vvol(ijk+jw)
         wscy3 = -c3y(ijk)/vvol(ijk+jw)
         wscz3 = -c3z(ijk)/vvol(ijk+jw)
!/vvol(ijk)
         wscx4 = -c4x(ijk)/vvol(ijk)
         wscy4 = -c4y(ijk)/vvol(ijk)
         wscz4 = -c4z(ijk)/vvol(ijk)
!/vvol(ijk)
         wscx5 = -c5x(ijk)/vvol(ijk+iw+kw)
         wscy5 = -c5y(ijk)/vvol(ijk+iw+kw)
         wscz5 = -c5z(ijk)/vvol(ijk+iw+kw)
!/vvol(ijk)
         wscx6 = -c6x(ijk)/vvol(ijk+iw+jw+kw)
         wscy6 = -c6y(ijk)/vvol(ijk+iw+jw+kw)
         wscz6 = -c6z(ijk)/vvol(ijk+iw+jw+kw)
!/vvol(ijk)
         wscx7 = -c7x(ijk)/vvol(ijk+jw+kw)
         wscy7 = -c7y(ijk)/vvol(ijk+jw+kw)
         wscz7 = -c7z(ijk)/vvol(ijk+jw+kw)
!/vvol(ijk)
         wscx8 = -c8x(ijk)/vvol(ijk+kw)
         wscy8 = -c8y(ijk)/vvol(ijk+kw)
         wscz8 = -c8z(ijk)/vvol(ijk+kw)
!
         do is = 1, nsp
!
            beta = 0.5*qom(is)*dt*transpos/clite
!
            cdotb = c1x(ijk)*bxv(ijk+iw) + (c1y(ijk)*byv(ijk+iw)+c1z(ijk)*bzv(&
               ijk+iw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+iw)**2+byv(ijk+iw)**2+bzv(ijk+&
               iw)**2))/vvol(ijk+iw)
!
            wscx1 = wscx1 - qom(is)*qdnv(ijk+iw,is)*(c1x(ijk)+beta*(c1y(ijk)*&
               bzv(ijk+iw)-c1z(ijk)*byv(ijk+iw)+beta*cdotb*bxv(ijk+iw)))*denom
!
            wscy1 = wscy1 - qom(is)*qdnv(ijk+iw,is)*(c1y(ijk)+beta*(c1z(ijk)*&
               bxv(ijk+iw)-c1x(ijk)*bzv(ijk+iw)+beta*cdotb*byv(ijk+iw)))*denom
!
            wscz1 = wscz1 - qom(is)*qdnv(ijk+iw,is)*(c1z(ijk)+beta*(c1x(ijk)*&
               byv(ijk+iw)-c1y(ijk)*bxv(ijk+iw)+beta*cdotb*bzv(ijk+iw)))*denom
!
            cdotb = c2x(ijk)*bxv(ijk+iw+jw) + (c2y(ijk)*byv(ijk+iw+jw)+c2z(ijk)&
               *bzv(ijk+iw+jw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+iw+jw)**2+byv(ijk+iw+jw)**2+&
               bzv(ijk+iw+jw)**2))/vvol(ijk+iw+jw)
!
            wscx2 = wscx2 - qom(is)*qdnv(ijk+iw+jw,is)*(c2x(ijk)+beta*(c2y(ijk)&
               *bzv(ijk+iw+jw)-c2z(ijk)*byv(ijk+iw+jw)+beta*cdotb*bxv(ijk+iw+jw&
               )))*denom
!
            wscy2 = wscy2 - qom(is)*qdnv(ijk+iw+jw,is)*(c2y(ijk)+beta*(c2z(ijk)&
               *bxv(ijk+iw+jw)-c2x(ijk)*bzv(ijk+iw+jw)+beta*cdotb*byv(ijk+iw+jw&
               )))*denom
!
            wscz2 = wscz2 - qom(is)*qdnv(ijk+iw+jw,is)*(c2z(ijk)+beta*(c2x(ijk)&
               *byv(ijk+iw+jw)-c2y(ijk)*bxv(ijk+iw+jw)+beta*cdotb*bzv(ijk+iw+jw&
               )))*denom
!
            cdotb = c3x(ijk)*bxv(ijk+jw) + (c3y(ijk)*byv(ijk+jw)+c3z(ijk)*bzv(&
               ijk+jw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+jw)**2+byv(ijk+jw)**2+bzv(ijk+&
               jw)**2))/vvol(ijk+jw)
!
            wscx3 = wscx3 - qom(is)*qdnv(ijk+jw,is)*(c3x(ijk)+beta*(c3y(ijk)*&
               bzv(ijk+jw)-c3z(ijk)*byv(ijk+jw)+beta*cdotb*bxv(ijk+jw)))*denom
!
            wscy3 = wscy3 - qom(is)*qdnv(ijk+jw,is)*(c3y(ijk)+beta*(c3z(ijk)*&
               bxv(ijk+jw)-c3x(ijk)*bzv(ijk+jw)+beta*cdotb*byv(ijk+jw)))*denom
!
            wscz3 = wscz3 - qom(is)*qdnv(ijk+jw,is)*(c3z(ijk)+beta*(c3x(ijk)*&
               byv(ijk+jw)-c3y(ijk)*bxv(ijk+jw)+beta*cdotb*bzv(ijk+jw)))*denom
!
            cdotb = c4x(ijk)*bxv(ijk) + (c4y(ijk)*byv(ijk)+c4z(ijk)*bzv(ijk))
            denom = dtcsq/(1. + beta**2*(bxv(ijk)**2+byv(ijk)**2+bzv(ijk)**2))/&
               vvol(ijk)
!
            wscx4 = wscx4 - qom(is)*qdnv(ijk,is)*(c4x(ijk)+beta*(c4y(ijk)*bzv(&
               ijk)-c4z(ijk)*byv(ijk)+beta*cdotb*bxv(ijk)))*denom
!
            wscy4 = wscy4 - qom(is)*qdnv(ijk,is)*(c4y(ijk)+beta*(c4z(ijk)*bxv(&
               ijk)-c4x(ijk)*bzv(ijk)+beta*cdotb*byv(ijk)))*denom
!
            wscz4 = wscz4 - qom(is)*qdnv(ijk,is)*(c4z(ijk)+beta*(c4x(ijk)*byv(&
               ijk)-c4y(ijk)*bxv(ijk)+beta*cdotb*bzv(ijk)))*denom
!
            cdotb = c5x(ijk)*bxv(ijk+iw+kw) + (c5y(ijk)*byv(ijk+iw+kw)+c5z(ijk)&
               *bzv(ijk+iw+kw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+iw+kw)**2+byv(ijk+iw+kw)**2+&
               bzv(ijk+iw+kw)**2))/vvol(ijk+iw+kw)
!
            wscx5 = wscx5 - qom(is)*qdnv(ijk+iw+kw,is)*(c5x(ijk)+beta*(c5y(ijk)&
               *bzv(ijk+iw+kw)-c5z(ijk)*byv(ijk+iw+kw)+beta*cdotb*bxv(ijk+iw+kw&
               )))*denom
!
            wscy5 = wscy5 - qom(is)*qdnv(ijk+iw+kw,is)*(c5y(ijk)+beta*(c5z(ijk)&
               *bxv(ijk+iw+kw)-c5x(ijk)*bzv(ijk+iw+kw)+beta*cdotb*byv(ijk+iw+kw&
               )))*denom
!
            wscz5 = wscz5 - qom(is)*qdnv(ijk+iw+kw,is)*(c5z(ijk)+beta*(c5x(ijk)&
               *byv(ijk+iw+kw)-c5y(ijk)*bxv(ijk+iw+kw)+beta*cdotb*bzv(ijk+iw+kw&
               )))*denom
!
            cdotb = c6x(ijk)*bxv(ijk+iw+jw+kw) + (c6y(ijk)*byv(ijk+iw+jw+kw)+&
               c6z(ijk)*bzv(ijk+iw+jw+kw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+iw+jw+kw)**2+byv(ijk+iw+jw+kw)&
               **2+bzv(ijk+iw+jw+kw)**2))/vvol(ijk+iw+jw+kw)
!
            wscx6 = wscx6 - qom(is)*qdnv(ijk+iw+jw+kw,is)*(c6x(ijk)+beta*(c6y(&
               ijk)*bzv(ijk+iw+jw+kw)-c6z(ijk)*byv(ijk+iw+jw+kw)+beta*cdotb*bxv&
               (ijk+iw+jw+kw)))*denom
!
            wscy6 = wscy6 - qom(is)*qdnv(ijk+iw+jw+kw,is)*(c6y(ijk)+beta*(c6z(&
               ijk)*bxv(ijk+iw+jw+kw)-c6x(ijk)*bzv(ijk+iw+jw+kw)+beta*cdotb*byv&
               (ijk+iw+jw+kw)))*denom
!
            wscz6 = wscz6 - qom(is)*qdnv(ijk+iw+jw+kw,is)*(c6z(ijk)+beta*(c6x(&
               ijk)*byv(ijk+iw+jw+kw)-c6y(ijk)*bxv(ijk+iw+jw+kw)+beta*cdotb*bzv&
               (ijk+iw+jw+kw)))*denom
!
            cdotb = c7x(ijk)*bxv(ijk+jw+kw) + (c7y(ijk)*byv(ijk+jw+kw)+c7z(ijk)&
               *bzv(ijk+jw+kw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+jw+kw)**2+byv(ijk+jw+kw)**2+&
               bzv(ijk+jw+kw)**2))/vvol(ijk+jw+kw)
!
            wscx7 = wscx7 - qom(is)*qdnv(ijk+jw+kw,is)*(c7x(ijk)+beta*(c7y(ijk)&
               *bzv(ijk+jw+kw)-c7z(ijk)*byv(ijk+jw+kw)+beta*cdotb*bxv(ijk+jw+kw&
               )))*denom
!
            wscy7 = wscy7 - qom(is)*qdnv(ijk+jw+kw,is)*(c7y(ijk)+beta*(c7z(ijk)&
               *bxv(ijk+jw+kw)-c7x(ijk)*bzv(ijk+jw+kw)+beta*cdotb*byv(ijk+jw+kw&
               )))*denom
!
            wscz7 = wscz7 - qom(is)*qdnv(ijk+jw+kw,is)*(c7z(ijk)+beta*(c7x(ijk)&
               *byv(ijk+jw+kw)-c7y(ijk)*bxv(ijk+jw+kw)+beta*cdotb*bzv(ijk+jw+kw&
               )))*denom
!
            cdotb = c8x(ijk)*bxv(ijk+kw) + (c8y(ijk)*byv(ijk+kw)+c8z(ijk)*bzv(&
               ijk+kw))
            denom = dtcsq/(1. + beta**2*(bxv(ijk+kw)**2+byv(ijk+kw)**2+bzv(ijk+&
               kw)**2))/vvol(ijk+kw)
!
            wscx8 = wscx8 - qom(is)*qdnv(ijk+kw,is)*(c8x(ijk)+beta*(c8y(ijk)*&
               bzv(ijk+kw)-c8z(ijk)*byv(ijk+kw)+beta*cdotb*bxv(ijk+kw)))*denom
!
            wscy8 = wscy8 - qom(is)*qdnv(ijk+kw,is)*(c8y(ijk)+beta*(c8z(ijk)*&
               bxv(ijk+kw)-c8x(ijk)*bzv(ijk+kw)+beta*cdotb*byv(ijk+kw)))*denom
!
            wscz8 = wscz8 - qom(is)*qdnv(ijk+kw,is)*(c8z(ijk)+beta*(c8x(ijk)*&
               byv(ijk+kw)-c8y(ijk)*bxv(ijk+kw)+beta*cdotb*bzv(ijk+kw)))*denom
!
         end do
!
!
         divwscx = c1x(ijk)*wscx1 + (c2x(ijk)*wscx2+(c3x(ijk)*wscx3+(c4x(ijk)*&
            wscx4+(c5x(ijk)*wscx5+(c6x(ijk)*wscx6+(c7x(ijk)*wscx7+c8x(ijk)*&
            wscx8))))))
!
         divwscy = c1y(ijk)*wscy1 + (c2y(ijk)*wscy2+(c3y(ijk)*wscy3+(c4y(ijk)*&
            wscy4+(c5y(ijk)*wscy5+(c6y(ijk)*wscy6+(c7y(ijk)*wscy7+c8y(ijk)*&
            wscy8))))))
!
         divwscz = c1z(ijk)*wscz1 + (c2z(ijk)*wscz2+(c3z(ijk)*wscz3+(c4z(ijk)*&
            wscz4+(c5z(ijk)*wscz5+(c6z(ijk)*wscz6+(c7z(ijk)*wscz7+c8z(ijk)*&
            wscz8))))))
!
         diag(ijk) = divwscx + divwscy + divwscz
!
      end do
!
      return
      end subroutine preconcr


      subroutine oldpreco(nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32&
         , a33, gx, gy, gz, gxtilde, gytilde, gztilde)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: a11(*)
      real(double) , intent(in) :: a12(*)
      real(double) , intent(in) :: a13(*)
      real(double) , intent(in) :: a21(*)
      real(double) , intent(in) :: a22(*)
      real(double) , intent(in) :: a23(*)
      real(double) , intent(in) :: a31(*)
      real(double) , intent(in) :: a32(*)
      real(double) , intent(in) :: a33(*)
      real(double) , intent(in) :: gx(*)
      real(double) , intent(in) :: gy(*)
      real(double) , intent(in) :: gz(*)
      real(double) , intent(out) :: gxtilde(*)
      real(double) , intent(out) :: gytilde(*)
      real(double) , intent(out) :: gztilde(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
!-----------------------------------------------
!
!
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
!
         gxtilde(ijk) = a11(ijk)*gx(ijk) + a12(ijk)*gy(ijk) + a13(ijk)*gz(ijk)
         gytilde(ijk) = a21(ijk)*gx(ijk) + a22(ijk)*gy(ijk) + a23(ijk)*gz(ijk)
         gztilde(ijk) = a31(ijk)*gx(ijk) + a32(ijk)*gy(ijk) + a33(ijk)*gz(ijk)
!
      end do
!
!
!
      return
      end subroutine oldpreco


      subroutine precond(nvtx, ijkvtx, a11, a12, a13, a21, a22, a23, a31, a32, &
         a33, gx, gy, gz, gxtilde, gytilde, gztilde)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer , intent(in) :: nvtx
      integer , intent(in) :: ijkvtx(*)
      real(double) , intent(in) :: a11(*)
      real(double) , intent(in) :: a12(*)
      real(double) , intent(in) :: a13(*)
      real(double) , intent(in) :: a21(*)
      real(double) , intent(in) :: a22(*)
      real(double) , intent(in) :: a23(*)
      real(double) , intent(in) :: a31(*)
      real(double) , intent(in) :: a32(*)
      real(double) , intent(in) :: a33(*)
      real(double) , intent(in) :: gx(*)
      real(double) , intent(in) :: gy(*)
      real(double) , intent(in) :: gz(*)
      real(double) , intent(out) :: gxtilde(*)
      real(double) , intent(inout) :: gytilde(*)
      real(double) , intent(inout) :: gztilde(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: r1, r2, r3
!-----------------------------------------------
!
!dir$ ivdep
!
      do n = 1, nvtx
!
         ijk = ijkvtx(n)
!
         r1 = gx(ijk)
         r2 = gy(ijk) - a21(ijk)*r1
         r3 = gz(ijk) - (a31(ijk)*r1+a32(ijk)*r2)
!
         gztilde(ijk) = r3*a33(ijk)
         gytilde(ijk) = r2*a22(ijk) - a23(ijk)*gztilde(ijk)
         gxtilde(ijk) = r1*a11(ijk) - (a12(ijk)*gytilde(ijk)+a13(ijk)*gztilde(&
            ijk))
!
      end do
!
!
      return
      end subroutine precond


      subroutine residue(ncells, ijkcell, nvtx, ijkvtx, iwid, jwid, kwid, ibp1&
         , jbp1, tiny, c1x, c2x, c3x, c4x, c5x, c6x, c7x, c8x, c1y, c2y, c3y, &
         c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z, c5z, c6z, c7z, c8z, vol, &
         phi, gradphix, gradphiy, gradphiz, vvol, dive, srce, rdt, scalsq, &
         residu)
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M

!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02
!...Switches: -yf -x1
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: ncells
      integer  :: nvtx
      integer  :: iwid
      integer  :: jwid
      integer  :: kwid
      integer  :: ibp1
      integer  :: jbp1
      real(double)  :: tiny
      real(double) , intent(in) :: rdt
      real(double) , intent(in) :: scalsq
      integer  :: ijkcell(*)
      integer  :: ijkvtx(*)
      real(double)  :: c1x(*)
      real(double)  :: c2x(*)
      real(double)  :: c3x(*)
      real(double)  :: c4x(*)
      real(double)  :: c5x(*)
      real(double)  :: c6x(*)
      real(double)  :: c7x(*)
      real(double)  :: c8x(*)
      real(double)  :: c1y(*)
      real(double)  :: c2y(*)
      real(double)  :: c3y(*)
      real(double)  :: c4y(*)
      real(double)  :: c5y(*)
      real(double)  :: c6y(*)
      real(double)  :: c7y(*)
      real(double)  :: c8y(*)
      real(double)  :: c1z(*)
      real(double)  :: c2z(*)
      real(double)  :: c3z(*)
      real(double)  :: c4z(*)
      real(double)  :: c5z(*)
      real(double)  :: c6z(*)
      real(double)  :: c7z(*)
      real(double)  :: c8z(*)
      real(double)  :: vol(*)
      real(double)  :: phi(*)
      real(double)  :: gradphix(*)
      real(double)  :: gradphiy(*)
      real(double)  :: gradphiz(*)
      real(double) , intent(in) :: vvol(*)
      real(double)  :: dive(*)
      real(double) , intent(in) :: srce(*)
      real(double) , intent(out) :: residu(*)
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: n, ijk
      real(double) :: rvvol
!-----------------------------------------------
!
!     calculate residue
!
      call gradf (nvtx, ijkvtx, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, c6x&
         , c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, c4z&
         , c5z, c6z, c7z, c8z, phi, gradphix, gradphiy, gradphiz)
!
      call axisgrad (ibp1, jbp1, iwid, jwid, kwid, gradphix, gradphiy, gradphiz&
         )
!
      do n = 1, nvtx
         ijk = ijkvtx(n)
         rvvol = 1./vvol(ijk)
         gradphix(ijk) = -gradphix(ijk)*rvvol
         gradphiy(ijk) = -gradphiy(ijk)*rvvol
         gradphiz(ijk) = -gradphiz(ijk)*rvvol
      end do
!
!
!     calculate div(grad(phi))
!
      call divc (ncells, ijkcell, iwid, jwid, kwid, c1x, c2x, c3x, c4x, c5x, &
         c6x, c7x, c8x, c1y, c2y, c3y, c4y, c5y, c6y, c7y, c8y, c1z, c2z, c3z, &
         c4z, c5z, c6z, c7z, c8z, vol, gradphix, gradphiy, gradphiz, dive)
!
!
!
      do n = 1, ncells
         ijk = ijkcell(n)
!
!
         residu(ijk) = (srce(ijk)-rdt*phi(ijk)-scalsq*dive(ijk))*vol(ijk)
!
      end do
!
      return
      end subroutine residue

	  subroutine vector_product(ax,ay,az,bx,by,bz,cx,cy,cz)
!
!		A routine to compute the vector product C = A x B
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      real(double) :: ax,ay,az,bx,by,bz,cx,cy,cz

	  cx=ay*bz-by*az
	  cy=-ax*bz+az*bx
	  cz=ax*by-ay*bx

	  end subroutine vector_product

	  	  subroutine vector_product_norm(ax,ay,az,bx,by,bz,cx,cy,cz)
!
!		A routine to compute the vector product C = (A x B)/B^2
!-----------------------------------------------
!   M o d u l e s
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double
      use modify_com_M
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      real(double) :: ax,ay,az,bx,by,bz,cx,cy,cz
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      real(double) :: bnorm
!-----------------------------------------------

	  bnorm=bx**2+by**2+bz**2
	  cx=(ay*bz-by*az)/bnorm
	  cy=(-ax*bz+az*bx)/bnorm
	  cz=(ax*by-ay*bx)/bnorm

	  end subroutine vector_product_norm
